class: center, middle, inverse, title-slide # Lecture 10 ## Multiple Groups ### Psych 10 C ### University of California, Irvine ### 04/20/2022 --- ## Hypotheses as models - When we have one independent variable that can take more than 2 values, comparing a Null model which assumes that every group has the same expected value with the Effects model which assumes that every group has a different expected value is not enough to answer our research question. -- - This is because the Effects model can be better than the Null in multiple ways. -- - For example, the Effects model will account for the data better if it is true that: -- 1. Every group is different `\(\mu_1 \neq \mu_2 \neq \mu_3\)`. 2. At least one group is different `\(\mu_1 = \mu_2 \neq \mu_3\)` (like in our students anxiety example). -- - Both of the previous statements are consistent with the main assumption of the Effects model, which in isolation can only tell us that at least one group is different. --- ## Hypothesis as models - As we showed with the anxiety example, we can reach an answer to our research question by comparing multiple models, which assume that one group is different from the rest. We do this with each group one at a time. -- - When our independent variable can take only 3 different values, this resulted in us having to look at the predictions and errors of 5 different models. -- - If we increase the number of groups (or values of our independent variable) this number will increase rapidly. For example, if we have 4 groups the number of models that we could compare is 9, and for 5 it grows to 17. -- - It is always better if we can reduce the number of models that we have to compare in order to answer a scientific question. -- - One way to achieve this goal is by making our research question more precise. --- ## A better approach to comparing anxiety levels - Going back to our anxiety levels example, instead of asking the question: -- - Are there any differences between the anxiety levels of First year students in the 218, 2019, and 2020 cohorts? -- - We could ask the question: -- - Are the anxiety levels of students from the 2020 cohort different from students in 2019 and 2018? -- - Notice that our new question is now centered around First year students on the 2020 cohort and not on differences between all cohorts. -- - In order to answer this question we only need to look at 3 models, first the Null, then the Effects and finally the 2020 model from last class. -- - Comparing Null and Effects model will tell us if there is at least one group that is different. Then our 2020 model will answer our actual question, Are the anxiety levels of First year students in the 2020 cohort different from students in the 2018 and 2019 cohorts? --- ## Hypothesis and models - Thus, just by redefining our experimental question we can reduce the number of comparisons that we have to make. -- - In the previous example, the second question makes more sense, as there is no real reason to believe that 2018 and 2019 cohorts should be different. However, due to the pandemic, we would expect First year students from the 2020 cohort to be more anxious (on average) than others. -- - One important thing to keep in mind is that we don't look for hypothesis that will reduce the number of models that we need to compare. A hypothesis should always come first. -- - In other words, we perform an experiment because we have a hypothesis regarding the outcomes of an experiment. This hypothesis will also inform our model comparisons. --- class: inverse, center, middle # Example ## New treatment for blood preasure --- # Research problem: - A pharmaceutical company is about to launch a new drug treatment for high blood pressure -- - Before that, they have to show that their treatment works better than a control substance (placebo) and is at least as good as the current standard of care (standard) treatment. -- - They have designed an experiment where participants are randomly assigned to one of three groups. There will be 25 participants assigned to each group for a total of 75. -- - Given our objective, which models would you consider for a comparison? .can-edit.key-likes[ Model 1: Model 2: Model 3: Model 4: Model 5: ] --- ## Experiment - Experimenters have recorded the blood pressure of each participant after 2 weeks of treatment, and obtained the following data:
--- ## Visualizing data - Another way in which we can visualize our data is by constructing a combination between a box-plot and a histogram: -- .pull-left[ ```r ggplot(data = drug_trial) + aes(x = condition) + aes(y = blood_pressure) + aes(color = condition) + geom_boxplot() + geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.7, method = "histodot",col = "#55555555",fill = "#55555555", binwidth = 2) + theme_classic() + xlab("Condition") + ylab("Blood pressure") + guides(fill = guide_legend("Exam"), color = "none") + theme(axis.title.x = element_text(size = 25), axis.title.y = element_text(size = 25)) ``` ] .pull-right[ <img src="data:image/png;base64,#lec-10_files/figure-html/hist-score-out-1.png" style="display: block; margin: auto;" /> ] --- ## Models - We will start by evaluating our two usual models, the Null and Effects model. -- - **Null model:** There are no differences in blood pressure levels between participants regardless of treatment condition. `$$y_{ij}\sim\text{Normal}(\mu,\sigma_0^2)$$` -- - **Effects model:** At least one of the three conditions has a different expected blood pressure level than the others. `$$y_{ij}\sim\text{Normal}(\mu_j,\sigma_e^2)$$` --- ## Null model - As in previous examples first we want to calculate the predictions of the model, the Sum of Squared Errors (SSE) and the Mean Squared Error of the Null. -- .pull-left[ ```r # sample size n_total <- nrow(drug_trial) # add null prediction and error drug_trial <- drug_trial %>% mutate("prediction_null" = mean(blood_pressure), "error_null" = (blood_pressure - prediction_null)^2) # calculate sum of squared error sse_null <- sum(drug_trial$error_null) # calculate mean squared error mse_null <- 1/n_total * sse_null # calculate bic bic_null <- n_total * log(mse_null) + 1 * log(n_total) ``` ] .pull-right[ - Prediction = 129.16 - SSE = 5135 - MSE = 68.47 - BIC = 321.29 ] --- ## Effects model - Now we have to calculate the predictions, SSE and MSE associated with the Effects model. Remember that the Effects model assumes that each group follows its own normal distribution: -- .pull-left[ ```r # predictions by group effects <- drug_trial %>% group_by(condition) %>% summarise("pred" = mean(blood_pressure)) # add null prediction and error drug_trial <- drug_trial %>% mutate("prediction_effects" = case_when(condition == "control" ~ effects$pred[1], condition == "standard" ~ effects$pred[3], condition == "new_drug" ~ effects$pred[2]), "error_effects" = (blood_pressure - prediction_effects)^2) # calculate sum of squared error sse_effects <- sum(drug_trial$error_effects) # calculate mean squared error mse_effects <- 1/n_total * sse_effects # calculate bic bic_effects <- n_total * log(mse_effects) + 3 * log(n_total) ``` ] .pull-right[ - Prediction: - control = 136.62 - standard = 130.03 - new drug = 120.83 - SSE = 1989.31 - MSE = 26.52 - `\(R^2\)` = 0.61 - BIC = 258.81 ] --- ## Treatment model - For this model we assume that the expectation of the treatment groups (not control) are different from the control group. We can express the model as: `$$y_{i2},y_{i3}\sim\text{Normal}(\mu,\sigma_2^2)$$` `$$y_{i1}\sim\text{Normal}(\mu_1,\sigma_2^2)$$` where `\(j=2\)` if the participant was in the **standard of care** condition and `\(j=3\)` if the participant was in the **new drug** condition. -- - This is the model that we are interested in! -- - As we have done with other models, we evaluate this one by calculating its predictions and error. --- ## Treatment model .pull-left[ ```r # create an indicator variable drug_trial <- drug_trial %>% mutate("id_treatment" = ifelse(test = condition == "control", yes = "control", no = "treatment")) # predictions by group treatment <- drug_trial %>% group_by(id_treatment) %>% summarise("pred" = mean(blood_pressure)) # add treatment prediction and error drug_trial <- drug_trial %>% mutate("prediction_treatment" = case_when(id_treatment == "control" ~ treatment$pred[1], id_treatment == "treatment" ~ treatment$pred[2]), "error_treatment" = (blood_pressure - prediction_treatment)^2) # calculate sum of squared error sse_treatment <- sum(drug_trial$error_treatment) # calculate mean squared error mse_treatment <- 1/n_total * sse_treatment # calculate bic bic_treatment <- n_total * log(mse_treatment) + 2 * log(n_total) ``` ] .pull-right[ - Prediction: - control = 136.62 - standard = 125.43 - new drug = 125.43 - SSE = 3047.31 - MSE = 40.63 - `\(R^2\)` = 0.41 - BIC = 286.47 ] --- ## Results - From our results we can see that the best model was the Effects, furthermore, the **treatment** model was also better at accounting for our observations than the Null model was. -- - From these comparisons we can safely answer the question that we had: Is the new drug better than a control substance and at least as good as the standard of care? -- - The comparison between the Null model and the Effects model shows that there is at least one group whose expectation is different from the rest. -- - From the comparison between treatment groups and control groups we know that treatments and the control condition are different (they have different means). -- - Finally, the comparison between the treatment and effects model shows us that it is better to assume that the 3 groups are different than assuming that the treatments are equally better than the control group. -- - Taken together, these results indicate that the 3 groups are different. --- ## Results - From our model comparisons we could conclude that: -- - Participants that used the new drug for 2 weeks had an average blood pressure of 120.83, while participants using the standard of care drug and a control substance had blood pressures of 130.03 and 136.62 respectively. This suggests that the new drug has an effect on the average blood pressure of people under treatment. -- - Given that the model that assumes that treatment conditions are different from the control is better than the Null and that the model that assumes that every group is different is better than the other two models, we can conclude that the new drug is better than the control group and better than the standard of care. -- - The Effects model accounted for 61% of the total variation in our observations. --- ## Visualizing model fit - Once we have selected a model we can visualize how close its predictions are to our observations. This is considered a visualization of the accuracy or fit of the model. -- - Although we can't use this as a method for model comparison, it gives the reader information regarding how good the model predictions are in the context of the experiment. -- - As we said at the beginning of the course, statistics are functions of random variables (our experimental outcomes) and as such they follow a probability distribution. -- - Our estimators for the expected value of each group (the average of each group) are statistics, and therefore, they follow their own probability distribution. -- - The good thing is that given the assumptions in our models (observations follow a Normal distribution) the distribution of our estimator is also Normal! --- ## Adding error bars to model predictions - To add error bars to our estimators (AKA confidence intervals) we have to take into account their distribution. -- - For the model that we selected using our model comparison approach, we will build an approximate confidence interval for prediction of each group using the following equation: `$$\hat{\mu}_j \pm 1.96 \times \sqrt{\frac{\hat{\sigma}^2}{n_j}}$$` -- - Where the value `\(1.96\)` comes from the normal distribution, and it's the value at which the normal accumulates `\(0.975\)` of its probability. `\(\hat{\mu}_j\)` is the prediction for group `\(j\)` and `\(\hat{\sigma}^2\)` is the Mean Squared eror of the best model. -- - The value of `\(n_j\)` is equal to the number of observations that we use to calculate each prediction. For example, with the Effects model we used 25 observations to calculate each `\(\hat{\mu}_j\)` which means that we would use `\(n_j = 25\)`. --- ## Adding error bars to model predictions - In comparison, with the new drug model we used only 2 predictions, so `\(n_{new}=25\)` for the new drug condition and `\(n_{other} = 50\)` for other conditions. Therefore, we need to calculate the sample size according to the model we selected. - In order to visualize our model's predictions in R first we have to create a new data with the variables that we need. -- ```r # create a new data frame with prediction by group and error confint_effects <- drug_trial %>% group_by(condition) %>% summarise("sample_size" = n(), "prediction" = mean(blood_pressure), "standard_error" = 1.96 * sqrt(mse_effects/sample_size)) ``` --- ## Visualizing model fit: .pull-left[ ```r ggplot(data = confint_effects) + aes(x = condition) + aes(y = prediction) + aes(color = condition) + geom_pointrange( aes(ymin = prediction - standard_error, ymax = prediction + standard_error, color = condition), position = position_dodge(0.3), fatten = 2, size = 2.8) + geom_dotplot(data = drug_trial, binaxis = 'y', stackdir='center', dotsize = 0.5, method = "histodot", col = "#555555", fill = "#555555", alpha = 0.3, binwidth = 2, mapping = aes(x = condition, y = blood_pressure)) + theme_classic() + xlab("Condition") + ylab("Blood pressure") + guides(fill = guide_legend("Exam"), color = "none") + theme(axis.title.x = element_text(size = 25), axis.title.y = element_text(size = 25)) ``` ] .pull-right[ <img src="data:image/png;base64,#lec-10_files/figure-html/eff-fit-out-1.png" style="display: block; margin: auto;" /> ] --- ## New drug model - If we want to do the same with the new drug model then we have to take into account that one of our predictions was obtained using 50 observations, which means that our estimate should be more precise (the more observations we have the more precise the estimator). ```r # create a new data frame with prediction by group and error confint_newdrug <- drug_trial %>% group_by(id_treatment) %>% summarise("sample_size" = n(), "prediction" = mean(blood_pressure), "standard_error" = 1.96 * sqrt(mse_treatment/sample_size)) ``` --- ## Visualizing model fit: New drug model .pull-left[ ```r ggplot(data = confint_newdrug) + aes(x = id_treatment) + aes(y = prediction) + aes(color = id_treatment) + geom_pointrange( aes(ymin = prediction - standard_error, ymax = prediction + standard_error, color = id_treatment), position = position_dodge(0.3), fatten = 2, size = 2.8) + geom_dotplot(data = drug_trial, binaxis = 'y', stackdir='center', dotsize = 0.5, method = "histodot", col = "#555555", fill = "#555555", alpha = 0.3, binwidth = 2, mapping = aes(x = id_treatment, y = blood_pressure)) + theme_classic() + xlab("Condition") + ylab("Blood pressure") + guides(fill = "none", color = "none") + theme(axis.title.x = element_text(size = 25), axis.title.y = element_text(size = 25)) ``` ] .pull-right[ <img src="data:image/png;base64,#lec-10_files/figure-html/new-fit-out-1.png" style="display: block; margin: auto;" /> ] --- ## Homework data link: ```r link <- "https://raw.githubusercontent.com/ManuelVU/psych-10c-data/main/homework3-p1.csv" ```